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A simple model which exhibits dynamical flame properties in ID is presented. It is investigated 
analytically and numerically. The results are applicable to problems of flame propagation in super- 
novae la. 



Observations of distant supernovae la explosions provide us with a very important information for modern cosmology 
([ll-Q)- probable scenario of type la supernovae is the following: the explosions are caused by a sufficiently fast 
thermonuclear burning of white dwarfs close to the Chandrasekhar limit (J^). After many years of atrophysical 
research, the accepted picture is that a supernova la explosion usually starts as a flame front, propagating slowly 
[ due to the thermal conduction (subsonic deflagration). Then the flame is transformed, due to a yet unknown reason, 
into a supersonic detonation wave (f^). The main observational properties of the supernovae la explosions must 
' depend significantly on the fiame propagation process, on the transition to detonation and on the instabilities which 
, are inherent to both slow and fast propagation regimes. These are presently far from being satisfactorily understood. 
They pose a challenge for theoretical physics, hydrodynamics and mechanics, because the supernovae la explosions 
cannot be investigated experimentally, and interiors of white dwarfs are not directly observable. 

In this paper we present and investigate, analytically and by means of numerical simulations, a very simple model 
of pulsating instability of a subsonic deflagration under conditions typical for white dwarfs. The instability occurs 

■ already in a one-dimensional geometry, so our model is one-dimensional. Because of the strong quantum degeneracy of 
the electron gas in white dwarfs, there is only a very weak expansion of the two-component electron-nuclei gas at the 
deflagration front. Combined with the one-dimensional geometry, this means that the motion of the fuel plays only a 

_ minor role, so we will consider the fuel as motionless. The reacting species in white dwarfs are nuclei, whereas heat is 
II transported by relativistic electrons and photons. This means that the Lewis number Le, characterizing the relative 
O role of the thermal conduction and fuel diffusion, is much higher than 1: Le 3> 1, and that the diffusion of reactants 
^ can be neglected. The only relevant dimensionless number, incorporated in our model, is the so-called Zeldovich 
^ number Ze. It characterizes the steepness of the reaction rates dependence on the temperature, or the thickness of the 
I I conductive region of the flame front relative to the effective thickness of a layer with reactions. Exact definitions of 
the Zeldovich number in the context of our model will be presented below. Usually the subsonic deflagration becomes 

■ pulsating for sufficiently high Ze. Our model may help to get insight into mechanisms of pulsations of the subsonic 
defiagration. It may also provide useful tests for more or less sophisticated numerical codes used for simulations of 
CO burning in white dwarfs. It cannot give, however, a final answer about the existence of pulsations in real flames. 

qI^ , Our main goal is to develop a fully analytically solvable model of burning in order to check whether the critical Ze 

■ number for the instability, predicted by the theory, coincides with that obtained in a numerical experiment. For the 
_ ^ ] sake of analytical solvability we will crudely simplify some properties of the matter in white dwarfs. Our model can 

I . be considered as a further simpliflcation of the model considered by . 

' Pulsating regimes of the subsonic deflagration have attracted considerable attention for a long time ( ) ■ Our 
I . attention to this problem was inspired additionally by a long-lived controversy between the results of |lCll | and [ll| . 

' In [ll| it was argued that the subsonic deflagration in white dwarfs undergoes a pulsating instability. This instabilit 
. ^ ] can have important consequences for the deflagration-to-detonation transition. However, numerical simulations of [1 

■ did not show any instability. The difference between these two results may be caused by different and complicated 
5h ' databases for nuclear reactions and/or by peculiarities of numerical methods. Our simple model, introduced below, 

. . . i may help separate these two quite different sources of the controversy. 

Here is a plan of the remainder of the paper. Sec. U is devoted to the formulation of our simple model. It also 
contains a traveling wave solution for arbitrary Ze. Linear stability of the traveling wave solution is investigated in 
Sec. ini We find there the critical Zeldovich number Zecr, below which the traveling wave is stable against small 
perturbations. Sec. IIIII is devoted to numerical simulations of the front propagation below and above the threshold. 
The results of this section make it clear that adequate codes are needed for simulations of the subsonic deflagration, 
and for a better understanding of nonlinear pulsating regimes of flame propagation. Our conclusions are presented in 
the last section of the paper. 
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I. THE MODEL 



Our model includes two dependent variables: the temperature T and the deflcient reagent fraction c. There are also 
two independent variables: time t and distance x. Neglecting diffusion of the deficient reagent (which is reasonable. 
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FIG. 1: Comparison of Arrhenius law for reaction rate and our step wise function for the same Ze = 9. 

because the kinetic coefficients in a white dwarf are such that Le 3> I), we introduce the foUowing system of equations 
for T and c: 



CdtT = ndlT + WujcQ{T - To), dtc = -ujcQ{T - To), 



(1) 



where 0(...) is the step-function. These equations describe deflagration burning in sohd propeUants, because two 
main physical processes which drive a slow front are present here: the thermal conductivity and the burning itself. 
The equations include the following constants: C is the thermal capacity of the fuel per unit volume; W is energy per 
unit volume of the deficient reagent fraction, released due to the reaction; and lj is the reaction rate, with dimensions 

Assuming that the deficient reagent fraction before ignition is equal to co, one can see that the fuel temperature 
after complete burnout will increase by the value 



Tf 



C 



(2) 



Assuming further that the temperature T(0) before ignition is much lower than Tf we will neglect the former one, 
setting T(0) = 0. 

In real white dwarfs, the thermal capacity depends on the temperature: for high density p ^ 10^ -f- 10^ g/cm^ all 
parameters are determined by the relativistic degenerate electrons, that is C oc T. But here we omit the dependence 
for the sake of possibility of analytical analysis. The latter concerns also the coefficient of thermal conduction k which 
is determined by both electrons and photons. 

The step wise approximation of the temperature dependence of nuclear reactions rates are used to have linear 
equations in both domains T < Tq and T > Tq. Intrinsic nonlinearity of the problem is moved only to boundary 
conditions between the two regions. This simplification will allow us to investigate the problem analytically. However 
we believe that our approximation of the temperature dependence represents to some extend very sharp temperature 
dependence of nuclear reaction rates, especially for high values of the Zeldovich number: Ze^ 1. See Fig. [T] 

Usually the Zeldovich number Ze is defined assuming Arrhenius law for the temperature dependence of reactions 
rates: w oc exp(— T^/T). In this case the Zeldovich number is defined usually as: Ze= Ta/Tj. To approximate this 
temperature dependence we may equate the temperatures for both dependence where the rates become e-fold lower 
than their maximum values at T = T/. Then we will have the following definition of Tq through Ta and Tf: 



Ze : 



Ta 

Tf 



To 



Tf~Tn 



or To = T/ 



T„ 



Tf+Ta 



= Tf 



Ze 



1 + Ze 



(3) 



Applying the same procedure to the power law for the reaction rate (w cx T"), we would obtain the following 
relationship: Tq = T/e-^/" « T/ (1 - n'^) (for n > 1). 

We may introduce now new dimensionless variables labelled by and defined as follows: 



t = Tt\ 



Ix : 



cqc; 



T : 



TfT, 



(4) 
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where 




K 



(5) 



Scales for x, t, T, c are chosen for all characteristic parameters of the traveling wave in dimensionless units to be 
equal to 1. This concerns the velocity of the traveling wave, the concentration c far before the wave, the temperature 
T far beyond the wave, and the characteristic width of the wave. We will use below only the dimensionless variables 
skipping the tilde 

Thus our dimensionless system is as follows: 

dtT^dlT + ujoceiT-To), dtc = -LOocQiT ^ To) , (6) 

where 

uo=iOT = To/{l-To); (0 < Tq < 1) . (7) 

Dimensionless value of the igniting temperature Tq is used in Eq. In initial physical units it is equal to T^/Tf. 
Our effective Zeldovich number [Eq. ([3])] in dimensionless units can be expressed as follows: 

Ze = ^o = T^. (8) 

The stationary traveling wave depends on x and t only via the combination — x — vt. For this reason it is 
convenient to introduce new spacial coordinates — x — vt instead of x, where v is the constant velocity of moving 
frame. Thus T and c depend now on t and ^. 

It is important that the system ([6]) is linear in both T > Tq and T < Tq. Nonlinearity appears only at the points, 
where the transition — >■ 1 in 0-function occur. Our equations ^ mean that T, dxT and c are continues at T = Tq. 
These conditions can be treated as matching conditions for solutions of the linear equations at T < Tq and T > Tq. 
As a result we may obtain the following system of equations instead of the system (jS]) for monotonic in ^ functions. 
It is convenient for the investigation of the stationary traveling wave and its linear stability: for > Tq 

dtT+ vdcT+ + 9|r+ + UJQC+, dtc+ = vd^c+ - ujqc+, (9) 

and for T_ < To 

dtT^ = vd^T^ + 9|T_, c_ = 1. (10) 

Here 

f(c f] - f f-(^^*) for e>e/W, 

^^^''^{f+itt) for C<^f{t), ^''> 

where / represents T or c, and ^/(i) is position of the front, where T = Tq. The matching conditions read as: 

T+{^f{t),t) = r_(e/(t),t)=To; 
d^T+{^f{t),t) = d^T.{^f{t),t); 

c+{^fit),t) = 1. (12) 

The stationary traveling wave obeys the following boundary conditions 

^^oo:T = 0, c=l; ^^-oo: dxT = 0, c = 0, (13) 

and the conditions: 

dtc = dtT^^f{t)=0. (14) 



The latter equality originates from a freedom due to the translation invariance of the problem, and hence from 
possibility to set the front at an arbitrary point of the ^-space. 



4 



The system ((9t- (|14p presents nonlinear eigen-value problem, with v being the eigen-value. There is unique solution 
of this problem: 

V = 1 

e>0: c=l, T = r_=Toe-«, 

C < : c = c+ = e"««, T = T+ = 1 — e"''^ . (15) 

LJQ + 1 

As a result we will set w = 1 as a constant velocity of the moving frame for the equations ^ and ([TUl) . 
We have for the traveling burning wave velocity in the initial physical units: 



c \n 



II. FLAME FRONT STABILITY 



We investigate here linear stability of the traveling wave solution obtained in the previous section. The latter one 
is called now as an unperturbed, and will be designated by the superscript T^he full solution is presented in the 

form 

r^r(o)+r(i)^ c = c(°)+cW, (16) 

where T^^) < T^o), c^^^ < c(o) and < 1. In this case our system can be linearized. Since the equations (|9]) 
and (|10p are linear from the very beginning, their linearization leads only to setting of the superscript "(1)", and to 
setting V = 1. However linearization of the matching conditions is not so trivial. They are transformed to: 

cJ^^a^Tf (o,t) + rf (o,i) = o, (17) 
^(i)5^yW(o^t) + ra)(o,t) = 0, (18) 

= S,fdlT^^\0,t)+d^T^}\0,t), (19) 
(0,<) + cW(0,i) = 0, (20) 

The boundary conditions p3l) give vanishing boundary conditions for t|_^'' and c^^ at ^ — > —oo and for at 
^ — +0O. There is no source for the perturbation of c before the front, so we set c_ = 0. 

Since Eqs. (O and ([TU)) and the matching conditions (|17p - (pn|) together with the boundary conditions at ^ — ±oo 
do not contain explicit dependence on time, their solution should have the following dependence on t: cx e^*, where 
p is an eigen-value of corresponding eigen-value linear boundary problem. Perturbation of the front position has the 
similar form. Thus 

OW-XeP*. (21) 

Eqs. ^ and PH)) do not contain even explicit dependence on i^. This property leads to a very simple form of the 
solutions at ^ > and ^ < 0: 

(22) 

P+"«)« , (23) 

^" ^^pt+(p+c^o)^ (24) 

(p + Uof + Wo 

where 

p=\^ + \, (25) 





— ae 




= I3e' 




= 7e' 



and necessary condition ReA < to obey the boundary conditions at +oo. 
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The presentation of solution (l?n) - (|24p contains 4 arbitrary constants x^ ct, /? and 7. However the solution should 
obey the matching conditions ([T7 | - (|20p . As a result we have homogeneous system of 4 linear equations for 4 variables. 
For the system of equations have nontrivial solution its determinant should be equal to 0. As a result we obtain 
the following equation for the eigen value p expressed trough A in accordance to Eq. (|25p under condition that 
{p + wqY + Wo 7^ 0: 

2 Wo A^ + {ujI + 2 Wo) A^ + 2 A + w^ 



(wo + 1) A2 + Wq + Wq 



0. (26) 



The cubic polynomial entering into the nominator of Eq. (|26p has a root corresponding to p = that could be 
assumed beforehand (see below) . As a result solution of Eq. is reduced to quadratic equation. Thus we have the 
following set of eigen values in terms of A: 



Ai 



Vw^~~8wJ^ + Wq 



A 



2 



^w^ - 8 Wo - wo 



A3 = -1. (27) 

The 3rd root A3 = —1 gives ps, = 0, a3 = — wo/(l + wo), = 1, 73 = 0. This eigen solution corresponds in according 
to Eq. (fT^ to the small shift of the stationary traveling front as a whole. Existence of such solution could be supposed 
in advance because of a translational symmetry of the problem. Thus physical roots are Ai_2- Real and imaginary 
parts of pi_2 versus wo are plotted in Figs. [2] and [3] If wo > 8, so that both Ai and A2 are real, then both eigen values 
Pi, 2 are positive. If, however, wo < 8, then 



Wo — 6wo ± i (wo — 2) ^ 8wo — Wg 

= ^ . 

Hence the traveling wave solution is stable against small perturbations at wq < 6 and unstable in the opposite case 
Wo > 6. The eigen values at the threshold wo = 6 are equal to pi. 2 = ±4i/\/3. Thus the perturbations at the threshold 
become purely oscillating. 

Expressing this result in terms of effective Zeldovich number Ze, introduced above, we may say that subsonic 
deflagration, having the form of the traveling wave, is stable in the frame of our model against small perturbations 
at Ze < Zecr — 6. There is no stable stationary moving flame front at Zc > Zccr = 6 in the frame of our model. 
Since time dependence of perturbations at the threshold Ze = Zecr are purely oscillating, one may assume that the 
subsonic deflagration at Ze > Zecr = 6 becomes oscillating. 
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FIG. 3: The complex growth rate p2 versus too = Ze— 1. 

III. NUMERICAL RESULTS 
A. The original model 

Such a simple system can be calculated numerically and it is a good test for analytical predictions. We use the 
Crank-Nicolson finite-difference method to solve the system 0. Initial and boundary conditions are set in two vifays. 
In the first a task of self-formation of the flame by a hot wall is set. In the whole region of calculation [0; L] cold 
unburncd matter is put, but the left wall is hot: 

T{t, X = 0) = 1, T{t, x = L) = 0, 
T(i = 0,a; e (0,L)) = 0, c(t = 0, x G (0, L)) = I. (28) 

If in this case a stationary front appears then it is natural for the system. The second way is to set distribution of 
the temperature and the concentration according to the analytical solution (jlSp . and to observe its evolution. 

Numerical methods suppose discritization of space dx and time dt. To obtain the physically correct solution 
we should properly choose these quantities. The first condition is dt = 10~^/a;o. This implies that only 1% of 
matter will burn on every numerical step, and therefore prohibits abrupt changes in the solution, and in such a way 
controls fluctuations. The second condition comes from the analytical solution. It follows from Eq. that there 
are two typical lengths in the system: 1 and I/wq. To resolve every change we need dx dxg = min(l,l/wo)- 
During numerical simulation several calculations were made with dx > dxo to determine a consequence of wrong 
discretization. 

The front coordinate x is determined by point c = 0.5 (the definition diverge with the theoretical definition T — Tq, 
but in the case of the stationary wave the points will be located on a constant distance from each other). The 
dependence of x{t) for the simulation with loq = 1 {dx <C dxo) is shown in Fig. 21 Good linear curve means that 
the front moves with a constant velocity and is stable. By fitting dx/dt the front velocity could be found and it is 
V = 1.00, what is a very good agreement with the theory. 

But it must be noted that despite good linear dependence, observed on large scales in Fig. HI small scales exhibit 
pulsations. These pulsations are fully a numerical effect, because their period is exactly T ~ dx/v (it can be easily 
tested by simulations with different dx). 

The evolution of front position x{t) when wq = 7 is shown in Fig. [5j According to the theory such regime should be 
unstable, what occurs in the simulation: a regime of moving with nonconstant velocity interchanges with the regime 
of the front standing. For the model ([5]), the front will not move after t = 20. This fact is interpreted below. 

Table |I] presents results for the first way of initial and boundary conditions set. The commentary column in 
the table shows propagation regime: "flame" means the fiame propagation, "therm" is the evolution of medium 
parameters like in Fig. [5] (which corresponds to the unstable regime). In this case the temperature undergo evolution 
like thermoconductivity without burning. Combinations of the terms "flame" and "term" in Table |T] mean that there 
are transition stages in these cases before establishing of an ultimate regime. The latter one is the last word in the 
combinations. Results for the analytical solution as an initial condition are shown in Table [TTl From those two tables 
we see that wq = 6 is a critical point for system parameters. When ujq < 6 a stationary front in a system can exist. 
When ujQ > 6 the front appears, but have nonconstant velocity and live for a limited period of time. The dependence 
x{t) for such front (in case of flrst boundary condition) is shown in Fig. [SJ 
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FIG. 4: The front position x{t) versus time t for cjo ~ 1- 
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FIG. 5: The front position x{t) versus time t for cjq ~ 7. 

Here we should emphasize simulations when dx > dxo is set. Two runs (wq = 1, dx = 1.5 and ujq = 4, = 0.5) 
show that a wrong dx results in a wrong behavior, such as an incorrect front velocity or an instability-like regime. 



TABLE I: Results for wall ignition: 





dx 


V 


comm. 


1.0 


0.05 


1.000 


flame 


4.0 


0.05 


0.996 


flame 


5.0 


0.05 


0.992 


flame 


5.5 


0.05 


0.993 


therm-flame 


5.8 


0.05 


0.993 


therm-flame 


6.0 


0.05 


6.15-M.24 


therm-flame 


6.5 


0.05 


5.56^0.99 


therm-flame-therm 


7.0 


0.05 


6.40-M.04 


therm-flame-therm 


8.0 


0.05 


4.32-1.08 


therm-flame-therm 


1.0 


1.5 


0.711 


flame 


4.0 


0.5 


4.74H-0.94 


therm-flame-therm 
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TABLE II: Analytic initial conditions: 



Wo 


dx 


V 


comm. 


1.0 


0.05 


1.000 


flame 


4.0 


0.05 


0.996 


flame 


0.0 


n no 
U.Uz 


i.UUD 


flame 


5.8 


0.02 


1.010 


flame 


6.0 


0.01 


1.019 


flame 


7.0 


0.01 




therm 


8.0 


0.01 




therm 


9.0 


0.01 




therm 


1.0 


1.5 


0.711 


flame 


4.0 


0.5 




therm 
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FIG. 6: Front proflles in terms of T versus x for wq = 6 (Tq ~ 0.857) at different time moments: to < t\ < t2 < ts. 

B. Modification of tlie simple model 

The model suffers from some nonpliysical effects: wlien T < Tq the burning rate drops to zero. This is why the 
flame stops at certain moment of time in unstable regime and never runs again (only heating by the left wall could 
ignite further burning). Lets consider a small modification of the burning rate: 

dtT = dlT + R{c,T), dtC^-R{c,T), (29) 

i?(c, T) = woce(r - To) + uj^c^QiTo - T)e(T), (30) 

with uji <C u)q. This modification allows burning at all temperatures, what is more physically correct. The condition 
LOi <C Wo imply that the model correction does not influence on the stationary flame and the previous theory. So when 
uj < uJo the flame spreads with a constant velocity. When ui > ujq the evolution changes: firstly according to previous 
calculations front decays and smoothes in the "thermoconductivity" regime, but after it a slow burning in region 
T < Tq raises the temperature to the critical value and the flame blaze up again. An example of such propagation for 
Wo = 8 is shown in Fig. [3 whereas the front position versus time is shown in Fig. [H 

This evolution is called front pulsations. Such a regime is a little bit different than "classical" pulsations: Vh — 
vq + wi sin at. Here the flame stops, blaze up and stops again. It moves by jerks. It is interesting that intervals 
between the jerks are considerably longer than the value |27r/Im p\ determined in the linear theory presented above. 
This relationship as well as the jerky motion of the front take place even for ujq that is only slightly higher than the 
threshold ujq — 6. 



9 



1.2 



1 

0.8 
H 0.6 
0.4 
0.2 


10 20 30 40 50 60 70 80 90 100 

X 

FIG. 7: Sequential front profiles in terms of T versus x at time moments: t\ < t2 < < t^. 
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FIG. 8: Front propagation in terms of x{t) for ujq = m modified model (solid line). For comparison dotted line displays 
moving with constant velocity v = 1. 

IV. DISCUSSION AND CONCLUSIONS 

We have presented here a simple analytically-solvable model for flame propagation which exhibits a different behavior 
depending on the dimensionless parameter Ze: the effective Zeldovich number expressed through our dimensionless 
parameter ljq. The theory is compared with numerical simulations, and a good agreement is observed. When Ze < 6 
both analytic and numerical solutions exhibit a constant-velocity front. Our analytical theory shows that the traveling 
flame front becomes unstable at Ze > 6, whereas our numerical simulations show that the front is destroyed and passes 
to jerk-like pulsations. 

We believe that, in the case of pulsating instability in a real flame, the characteristic behavior will be the same as in 
this model, but the critical Ze can be different. Therefore, the proposed model can mimic the behavior of numerical 
methods, and serve as a test for them. Our present work cannot say anything about the existence of such pulsations 
in real white dwarfs. However, if such a pulsating jerk-like regime of slow flame propagation indeed takes place in 
white dwarfs, then it could be able to trigger the transition to detonation. 
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